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ABSTRACT 

To explore the connection between the global physical properties of galaxies and their far-infrared 
(FIR) spectral energy distributions (SEDs), we study the variation in the FIR SEDs of a set of 
hydrodynamically simulated galaxies that are generated by performing dust radiative transfer in post¬ 
processing. Our sample includes both isolated and merging systems at various stages of the merging 
process and covers infrared (IR) luminosities and dust masses that are representative of both low- and 
high-redshift galaxies. We study the FIR SEDs using principle component analysis (PCA) and find 
that 97% of the variance in the sample can be explained by two principle components (PCs). The first 
PC characterizes the wavelength of the peak of the FIR SED, and the second encodes the breadth of 
the SED. We find that the coefficients of both PCs can be predicted well using a double power law in 
terms of the IR luminosity and dust mass, which suggests that these two physical properties are the 
primary determinants of galaxies’ FIR SED shapes. Incorporating galaxy sizes does not significantly 
improve our ability to predict the FIR SEDs. Our results suggest that the observed redshift evolution 
in the effective dust temperature at fixed IR luminosity is not driven by geometry: the SEDs of 
z ~ 2 — 3 ultraluminous IR galaxies (ULIRGs) are cooler than those of local ULIRGs not because the 


high-redshift galaxies are more extended but rather 
luminosity. Finally, based on our simulations, we 
that depend on both IR luminosity and dust mass. 

1. INTRODUCTION 

Understanding what drives the variation in the far- 
infrared (FIR) spectral energy distributions (SEDs) of 
galaxies is a key goal if we wish to maximize the physical 
insight that can be gleaned from the wealth of data that 
is rapidly being collected in the era of Herschel (Pilbratt 
et al. 2010) and the Atacama Large Millimeter Array 
(ALMA). Unfortunately, interpreting FIR SEDs is dif¬ 
ficult because spatially resolved data are not available 
for the vast majority of galaxies (although ALMA will 
help greatly in this regard). Moreover, there are various 
degeneracies (e.g., between the dust temperature distri¬ 
bution and the power-law index of the dust emissivity 
curve; see Kelly et al. 2012 and references therein) that 
are difficult to break. 

FIR SEDs are fit using a wide variety of methods (see 
Casey et al. 2014 for a recent review). Often, one or 
a sum of a few modified blackbody SEDs are used to 
fit FIR SEDs, but the physical information that can be 
gained from such models is limited (see, e.g., Hayward 
et al. 2012 and Smith et al. 2013 for detailed discussions). 
Another common approach is to use templates: by fitting 
the available photometry with an IR SED template, the 
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because they have higher dust masses at fixed IR 
introduce a two-parameter set of SED templates 


IR luminosity can be estimated. Sometimes one or a 
few SEDs of well-studied IR luminous galaxies, such as 
Arp 220, are used. Various empirically based, single¬ 
parameter templates are also widely used (e.g.. Chary & 
Elbaz 2001; Pope et al. 2008; Rieke et al. 2009; Magdis 
et al. 2012). 

Although empirically based FIR SED templates have 
many practical uses, physical models are necessary if one 
wishes to learn about, e.g., the radiation field and dust 
properties of a galaxy. Because of the computational ex¬ 
pense of radiative transfer and difficulty of constraining 
many free parameters with a limited number of FIR pho¬ 
tometric data points, approaches to modeling FIR SEDs 
without doing radiative transfer calculations have been 
advocated (e.g.. Desert et al. 1990; Devriendt et al. 1999; 
Dale et al. 2001; Dale & Helou 2002; Draine & Li 2007; 
da Cunha et al. 2008; Somerville et al. 2012). Some au¬ 
thors (e.g.. Dale et al. 2001; Dale & Helou 2002) model 
FIR SEDs by assuming that the mass distribution of dust 
that is exposed to radiation fields of different intensities 
can be described by a truncated power law. Draine & 
Li (2007) expanded this model by adding a delta func¬ 
tion to the intensity distribution (at the minimum inten¬ 
sity) to represent diffuse ISM dust (see also Draine et al. 
2007). Such models can provide good fits to observed 
SEDs. Similarly, other authors (e.g., Kovacs et al. 2010; 
Magnelli et al. 2012) parameterize SEDs by assuming a 
power-law distribution of dust temperatures. 

Radiative transfer calculations, either with simplified 
assumptions about the geometry of dust with respect to 
stars (Witt & Gordon 1996, 2000; Silva et al. 1998; Ef- 
stathiou et al. 2000; Gordon et al. 2001; Takagi et al. 
2003; De Geyter et al. 2014) or active galactic nuclei 
(AGN; e.g., Fritz et al. 2006; Siebenmorgen & Krugel 
2007; Stalevski et al. 2012), or more complicated geome¬ 
tries (Dopita et al. 2005; Popescu et al. 2011; De Looze 
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et al. 2012, 2014; De Geyter et al. 2015), have been qual¬ 
itatively successful in producing SEDs similar to those 
observed. It is also possible and (has become increas¬ 
ingly common) to ‘forward-model’ FIR SEDs by per¬ 
forming three-dimensional (3D) dust radiative transfer 
on the outputs of hydrodynamical simulations of galax¬ 
ies in post-processing (e.g. Jonsson 2006; Jonsson et al. 
2006, 2010; Chakrabarti et al. 2007, 2008; Chakrabarti 
& Whitney 2009; Narayanan et al. 2010a, b; Hayward 
et al. 2011, 2012, 2014b; Snyder et al. 2013; Dominguez- 
Tenreiro et al. 2014; Lanz et al. 2014; Martmez-Galarza 
et al. 2014; Granato et al. 2015). 

Because explicit radiative transfer calculations can 
best capture complicated source and dust geometries; di¬ 
rectly treat dust absorption (including not only absorp¬ 
tion of primary radiation from stars and AGN but also 
dust self-absorption), scattering, and re-emission; and 
fully characterize the 3D distribution of dust temper¬ 
ature, which depend on both the local radiation field 
and the grain properties, they provide the best tool for 
studying how FIR SEDs depend on galaxy properties. 
By performing 3-D Monte Carlo radiative transfer calcu¬ 
lations for idealized geometries (i.e., not outputs of hy¬ 
drodynamical simulations), Misselt et al. (2001) studied 
the dependence of the shape of the FIR SED on global 
parameters of the emitting regions. They found that 
higher dust mass leads to colder SEDs when the other 
parameters in their model are fixed. This is a simple 
consequence of thermal equilibrium. Moreover, a ‘shell’ 
(aka foreground screen) geometry results in a broader 
SEDs compared to a geometry in which stars and dust 
are mixed (which they term the ‘dusty’ geometry) be¬ 
cause of differences in the temperature structures of the 
two models. 

Performing radiative transfer on hydrodynamical sim¬ 
ulations of galaxy formation has the benefit of including 
more realistic source and dust geometries than simpler 
approaches, such as that of Misselt et al. (2001). This 
approach has been demonstrated (Jonsson et al. 2010) to 
yield SEDs that agree well with observed local samples 
of normal star forming galaxies from the Spitzer Infrared 
Nearby Galaxies Survey (SINGS; Kennicutt et al. 2003) 
and local luminous infrared (IR) galaxies (LIRGs) from 
the Great Observatories All-sky LIRG Survey (GOALS; 
Armus et al. 2009). Moreover, Lanz et al. (2014) have 
shown that the SEDs of local interacting galaxies from 
the Spitzer Interacting Galaxies Survey (SIGS; Lanz 
et al. 2013; Brassington et al. 2015) can be fit reasonably 
well with SEDs predicted in this manner; this is also the 
case for 24-/rm-selected galaxies at z ~ 0.3—2.8 (Roebuck 
et ah, in preparation). Unfortunately, performing 3D ra¬ 
diative transfer on hydrodynamical simulations is orders- 
of-magnitude more computationally expensive than less 
detailed approaches that assume smooth axisymmetric 
geometries (e.g., Silva et al. 1998): the former typically 
requires of order 10 — 10^ CPU-hours per galaxy, whereas 
the latter requires at most a few minutes. Thus, it is in¬ 
structive to determine whether the high computational 
expense of performing 3D dust radiative transfer on hy¬ 
drodynamical simulations of galaxies can be avoided, as 
it may be possible to characterize the variation among 
the simulated SEDs using only a few global parameters. 

In this work, we investigate whether we can predict the 


FIR SEDs of simulated galaxies, which were calculated 
using 3D dust radiative transfer in previous works, us¬ 
ing a simple, computationally inexpensive method. The 
method that we use to study these SEDs is principle com¬ 
ponent analysis (PGA). In our case, PGA yields FIR SED 
eigenvectors (principle components, hereafter PCs) such 
that linearly combining them with the mean SED of the 
sample can capture the variance in the simulated SEDs. 
The coefficients of each PC are different for each galaxy, 
and studying how the coefficients depend on global pa¬ 
rameters, such as the star formation rate (SFR) or IR 
luminosity, can give physical intuition regarding what 
drives the dispersion in the FIR SEDs. 

The simulated galaxy SEDs used in this work have 
been analyzed extensively in previous works. These or 
similar simulations have been demonstrated to exhibit 
good agreement with the SEDs/colors of various classes 
of real galaxies, including local ‘normal’ galaxies (Jons¬ 
son et al. 2010) and (U)LIRGs (Younger et al. 2009; Jon¬ 
sson et al. 2010; Lanz et al. 2014), high-redshift dusty 
star-forming galaxies (Narayanan et al. 2010a,b; Hay¬ 
ward et al. 2011, 2012; Roebuck et ah, in prep.), obscured 
AGN (Snyder et al. 2013; Roebuck et ah, in preparation), 
post-starburst galaxies (Snyder et al. 2011), and compact 
quiescent galaxies (Wuyts et al. 2010). Thus, although 
the simulations used herein naturally have some associ¬ 
ated limitations (see Section 8.6), they have the advan¬ 
tage of being well-tested and in agreement with many 
observational constraints. Furthermore, to the best of 
our knowledge, they still represent the state of the art in 
terms of ultraviolet (UV)-millimeter (mm) SEDs gener¬ 
ated from 3D hydrodynamical simulations of galaxies. 

In addition to determining what physical insights 
about galaxies can be gained from their FIR SEDs, we 
hope to identify possibilities for improving how semi- 
analytic models (SAMs) of galaxy formation predict FIR 
SEDs. Some SAMs (e.g., Silva et al. 1998; Granato 
et al. 2000) employ radiative transfer calculations that 
assume axisymmetric geometries or analytic models that 
are designed to capture the results of such calculations 
(e.g., Gonzalez et al. 2011). Such calculations have been 
widely employed, but it is unclear how well their results 
agree with those of 3D radiative transfer calculations per¬ 
formed on hydrodynamical simulations of galaxies with 
more complex geometries. Other SAMs (e.g., Somerville 
et al. 2012) use empirically derived templates that are a 
function of a single parameter, Lir, to predict FIR SEDs. 
However, such an approach is insufficient; for example, it 
cannot capture the redshift evolution of the effective dust 
temperature~IR luminosity relation (Casey et al. 2014 
and references therein) by construction. Thus, we aim 
to determine one or more additional physical parameters, 
besides Lir, that can be used to predict the FIR SEDs 
of galaxies. Having determined what additional param- 
eter(s) is necessary to predict the variation in IR SEDs, 
we can then generate a set of multi-parameter SED tem¬ 
plates. By incorporating these templates in a SAM, we 
may be able to resolve the tension between the observed 
FIR and submillimeter number counts and those pre¬ 
dicted by the SAM (Somerville et al. 2012; Niemi et al. 
2012); this will be explored in a future work. These tem¬ 
plates could be also used in semi-empirical models (e.g., 
Bethermin et al. 2012, 2013; Bernhard et al. 2014). Fi¬ 
nally, these two-parameter templates will be useful for 
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fitting observed SEDs to infer the total IR luminosity 
and predict fluxes at wavelengths for which data are not 
available. 

The remainder of this paper is organized as follows: 
in Section 2, we describe the properties of the simulated 
galaxy SED dataset that we use and summarize the de¬ 
tails of the hydrodynamical simulation and dust radia¬ 
tive transfer calculation methods. Section 3 summarizes 
the PCA technique and how we use the PCA results to 
predict the simulated galaxies’ SEDs. In Section 4, we 
present the results of the PCA of the simulated galaxy 
SEDs. Section 5 demonstrates that the dust mass is a 
key parameter, in addition to the IR luminosity, that de¬ 
termines the SED shape, whereas in Section 6, we show 
that incorporating the galaxy size does not improve our 
ability to predict the EIR SEDs. In Section 7, we present 
a two-parameter family of SED templates that depend 
on the IR luminosity and dust mass. In Section 8, we 
discuss some observational support for the importance 
of dust mass in determining the SED shape; speculate 
regarding the unimportance of galaxy size, the origin of 
catastrophic failures to predict the SEDs, and the pos¬ 
sible implications of using the proposed two-parameter 
templates in SAMs; and highlight some limitations of 
our work. Section 9 lists our primary conclusions. 

2. SIMULATED GALAXY SED DATASET 

In this work, we analyze two sets of EIR SEDs of 
simulated isolated and merging disk galaxies that were 
presented in previous works. We will first summarize 
the properties of the two datasets and then briefly dis¬ 
cuss the simulation methodology. The first SED dataset 
was originally presented in Lanz et al. (2014) and Hay¬ 
ward et al. (2014b). The four progenitor galaxies span a 
stellar (baryonic) mass range of 6 x 10® — 4 x 10^° Mq 
(10® —5x 10^® Mq), and their properties were selected to 
be representative of typical star-forming galaxies in the 
z ~ 0 universe (see Cox et al. 2008 for details). Each 
of the four progenitors was simulated in isolation, and 
binary mergers of all possible combinations of progeni¬ 
tors were simulated for a single representative orbit (i.e., 
the results for other ‘non-special’ orbits are similar; the 
results for perfectly co-planar and other rare special con¬ 
figurations can sometimes differ significantly; e.g., Cox 
et al. 2008), thereby yielding 10 merger simulations. The 
total dataset contains ^ 12, 000 SEDs. We refer to this 
dataset as the z ~ 0 dataset. 

The second set of SEDs of simulated isolated disk and 
merging galaxies was presented in Hayward et al. (2011, 
2012, 2013b). For this dataset, the structural properties 
of the progenitor (bulgeless) disk galaxies were scaled to 
z ~ 3 following the method of Robertson et al. (2006), 
and the initial gas fractions of the disks, 0.6-0.8, are sig- 
nihcantly greater than those of the z ~ 0 simulations.^ 
Because the original purpose of these simulations was to 
model high-redshift submillimeter galaxies (SMGs), the 
progenitor disks span a relatively narrow baryonic mass 
range of^ 1 — 4x10^^ Mq, but a variety of merger orbits 
and mass ratios are included (see Hayward et al. 2012 for 

^ Such large gas fractions were used to ensure that the gas frac¬ 
tion at the time of coalescence was consistent with observational 
constraints for z ~ 2 galaxies. For all snapshots considered in this 
work, the gas fractions are < 0.5. See Hayward et al. (2013b) for 
more details. 


details). This dataset contains 37 hydrodynamical sim¬ 
ulations, from which ~ 6500 SEDs were calculated. We 
refer to this dataset as the z ^ 3 dataset. 

The SEDs contained in the two datasets were gen¬ 
erated by performing dust radiative transfer in post¬ 
processing on the outputs of 3-D hydrodynamical sim¬ 
ulations. The full methodology is described in the afore¬ 
mentioned works and references therein, so we will only 
summarize it here. Eirst, idealized isolated (i.e., non- 
cosmological) galaxy models were created following the 
method described in Springel et al. (2005). Each initial 
disk galaxy is composed of a dark matter halo, stellar and 
gaseous disks, and a supermassive black hole (SMBH); 
for the z ~ 0 simulations only, a stellar bulge is also in¬ 
cluded. Then, the isolated galaxies and binary mergers 
of these galaxies were simulated using a heavily modi¬ 
fied version of the gadget- 2 Wbody/smoothed-particle 
hydrodynamics (SPH) code (Springel 2005).® The sim¬ 
ulations inclnde the effects of gravity, hydrodynamical 
interactions, and radiative heating and cooling. Star for¬ 
mation and stellar feedback are incorporated via the two- 
phase sub-resolution interstellar medium (ISM) model of 
Springel & Hernquist (2003), and BH accretion and AGN 
feedback are treated following Springel et al. (2005). 
Each gas particle is enriched with metals according to 
its associated SER, assuming a yield of 0.02. 

Every 10 Myr, ‘snapshots’ of the physical state of 
the simulations were saved. Then, to calculate UV- 
mm SEDs, 3-D dust radiative transfer was performed 
in post-processing on a subset of the snapshots using the 
code SUNRISE (Jonsson 2006; Jonsson et al. 2010). Eor 
a given snapshot, the SUNRISE calculation proceeds as 
follows: the stellar and BH particles in the GADGET-2 
simulation, which are the sources of radiation, are as¬ 
signed source SEDs according to their properties (age 
and metallicity for the star particles and luminosity for 
the BH particles). The metal distribution from the simu¬ 
lation is projected onto an octree grid for the purpose of 
calculating the dust optical depths (assuming a dust-to- 
metal density ratio of 0.4; e.g., Dwek 1998; James et al. 
2002; Sparre et al. 2014). The Milky Way Ry = 3.1 dust 
model of Draine & Li (2007) was used. 

The most significant uncertainty (of which we are 
aware) in these calculations is the sub-resolution struc¬ 
ture of the ISM. Specifically, the simulations do not re¬ 
solve the ISM structure on scales of less than a few hun¬ 
dred parsecs, but observations clearly indicate that the 
ISM has substantial structure on smaller scales. We can 
crudely investigate this uncertainty through the use of 
two alternate treatments for the sub-resolution ISM.® In 
the first treatment (referred to as ‘multiphase on’ or ‘de¬ 
fault ISM’ in previous works), it is assumed that the 
cold clouds implicit in the Springel & Hernquist (2003) 
model have negligible volume filling factors. Thus, for 
the purpose of the radiative transfer calculations, this 

® The traditional density-entropy formulation of SPH was used. 
However, we note that the results of idealized non-cosmological 
galaxy simulations such as these are relatively insensitive to the 
numerical method (Hayward et al. 2014a) and thus do not consider 
the use of traditional SPH to be a significant limitation. 

® This issue has been discussed extensively in previous works 
(e.g., Hayward et al. 2011; Snyder et al. 2013; Lanz et al. 2014; 
Hayward &; Smith 2015), and we refer the interested reader to 
those papers for additional details. 
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dust is ignored. To account for the obscuration of young 
stars from their birth clouds, a sub-resolution model for 
HII and photodissociation regions (Groves et al. 2008) is 
used. In the second treatment (‘multiphase off’ or ‘al¬ 
ternate ISM’), the entire dust mass (i.e., dust in both 
the implicit diffuse phase and cold clouds) is considered 
in the radiative transfer calculations. In this work, we 
use the ‘multiphase off’ model to ensure that all of the 
dust emission is calculated self-consistently (i.e., there 
is no component from the sub-resolution HII and pho¬ 
todissociation region model). However, we have checked 
that our conclusions are qualitatively unchanged if the 
‘multiphase on’ model is used. 

After the source and dust properties are specified, 
radiation transfer is performed using the Monte Carlo 
method to calculate the effects of dust absorption, 
scattering, and re-emission. Importantly, dust self¬ 
absorption is treated using an iterative procedure in 
which the transfer of IR radiation and dust temperature 
calculation are repeated until the luminosity absorbed 
in each cell is sufficiently converged. For each snapshot, 
this process yields spatially resolved UV-mm SEDs of the 
simulated galaxy/merger viewed from 7 viewing angles. 

We only analyze the SEDs of galaxies that are either 
isolated systems or mergers that are experiencing coales¬ 
cence (defined by a black hole separation of <1 kpc) or 
are post-coalescence. The reason for this choice is that 
we only want to study systems that can be considered 
a single galaxy. Using this black hole separation crite¬ 
rion excludes early stage mergers, which would often be 
considered separate systems in low-redshift observations 
(although they would be unresolved in the FIR at high 
redshift when observed with single-dish FIR or submil¬ 
limeter telescopes). Moreover, when the system consists 
of widely separated galaxies, the radiative transfer within 
individual galaxies is decoupled (i.e., the radiation within 
one galaxy does not have a significant effect on the dust 
in the other galaxy when the two galaxies are separated 
by many kiloparsecs). Using a less conservative criterion, 
such as 5 kpc, does not qualitatively affect our results. 


3. PREDICTING THE SEDS BASED ON PCA 

Given a sample of FIR SEDs, PCA finds the mean SED 
of the entire sample ((AL^)) and based on the mean SED, 
finds the PCs that encapsulate most of the variance of the 
entire sample. PCs are eigenvectors that are orthogonal 
to each other. The PCs are not guaranteed to have a 
physical interpretation; rather, they are a tool to reduce 
the dimensionality of the problem. The PCs are ranked 
in terms of how much of the variance is explained by each 
PC. 

Because we are interested in studying the shape of the 
SEDs, we normalize the EIR SEDs by dividing them by 
the IR luminosity. The units of the data on which the 
PCA is performed do not affect the results, and we have 
chosen to work in dimensionless units of ALa/Tir- 

The SED belonging to a given galaxy (j) in the sample 
can be reconstructed as a linear combination of principle 


We define Lir as the integral of the SED over the wavelength 
range of 8-1000 fim. 



Figure 1. The black solid line corresponds to the mean FIR SED 
of our full sample; it has been normalized by dividing by Lir. The 
first (PCI; blue dashed line) and second (PC2; green dash-dotted 
line) PCs are also shown. These PCs are added to the mean with a 
unique coefficient for each galaxy to reconstruct the individual FIR 
SEDs. The fraction of the variance in the data that is explained 
by each PC is indicated as a percentage in the legend. Together, 
the top 2 PCs can explain 97% of the total variance of the data in 
the full SED sample, which includes both low- and high-redshift 
simulated galaxies. 


components: 

N 

XL\ j = (ALa) + 'y ] Cjj X T’C'i, (I) 

i=l 

where N is the number of PCs used, PCi is the ith PC, 
and Cij is the coefficient of the Ah PC for the jth galaxy. 

Having identified the most important PCs, we then 
examine how the coefficients of the PCs correlate with 
various global physical parameters of the galaxies. Our 
goal is to be able to predict the PC coefficients, and thus 
the FIR SED, of a galaxy based on a small number of 
galaxy properties, such as the SFR. To predict the PC 
coefficients, we use functions of the following form: 

M 

C'ij = Pi,k log Pt,j,k , (2) 

k=l 

where C'^ is the predicted value of the ith coefficient for 
galaxy j, ai is the fit intercept for the ith PC coefficient, 
Pi,k is the fit coefficient for the *th PC coefficient and 
fcth parameter, and Pi,j,k is galaxy j’s value for the fcth 
physical property used to predict the *th PC coefficient; 
example properties include the SFR and dust mass. As 
we will discuss below, we found that 2 PCs were sufficient 
to explain 97% of the variance. Thus, we used N = 2. 
To predict the PC coefficients, we tried relations with 
both a single physical property {M = 1) and a pair of 
physical properties (M = 2). 

Given predicted values for a galaxy’s PC coefficients 
determined using Equation (2), we predict its SED as 
follows: 

N 

AL'a,, = (ATa)+E^1. (3) 

To quantify how well an SED can be predicted, we use 
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Figure 2. The top set of four figures shows how the coefficient 
of the first PC depends on four global physical parameters of the 
galaxies in the entire sample: IR luminosity {upper left), AGN 
luminosity {upper right), SFR {lower left), and dust mass {lower 
right). The bottom set of four figures shows the same plots for C2. 
In each panel, the parameter value is shown on the y-axis, and the 
ic-axis corresponds to the value of the coefficient. The color of each 
bin corresponds to the logarithm of the number of points in the bin, 
as specified by the colorbars. Cl is correlated with IR luminosity, 
SFR and AGN luminosity, although there is a large scatter. This 
result indicates that as the IR luminosity, SFR or AGN luminosity 
is increased, the SED peak shifts toward shorter wavelengths. G2 
is anti-correlated with all of the four physical parameters, although 
the correlations are weaker (see the text for details). This result in¬ 
dicates that there is a weak tendency for higher-luminosity or/and 
higher-dust mass galaxies to have narrower FIR SEDs. 


the reduced chi squared value: 


X 


2 

r 


1 

df h 


( 4 ) 


where where L\^ and denote the true and predicted 
luminosity density values at wavelength Ai, is the 
uncertainty at wavelength Ai, P = 20 is the total number 
of wavelength bins in the FIR, and df = P — M — 1 is 
the number of degrees of freedom. Because our data are 
noise-free, for the purposes of calculating the reduced 
value, we have arbitrarily assumed a signal-to-noise 
ratio of 5 in each band (i.e., a\^ = 0.2Lx.). Thus, the 


Xr values are useful in a relative sense (i.e., to determine 
which SEDs are predicted better than others), but the 
absolute values are not very meaningful. 

4. PCA RESULTS 

We performed PCA on the FIR (A > 25 /rm) SEDs 
from our entire sample, which includes both low- and 
high-redshift simulated galaxies. Figure 1 shows the 
mean SED and first two PCs of our sample. The percent¬ 
age of the variance in the SEDs that can be explained by 
each PC is indicated in the figure. Each PC is multi¬ 
plied by a coefficient (which is unique for each galaxy) 
and added to the mean SED of the sample to reconstruct 
the FIR SED of a particular galaxy. The coefficient can 
be negative or positive. From Figure 1, it is clear that 
adding the first PC (PCI) with a positive coefficient to 
the mean FIR SED tends to make the SED warmer with 
respect to the mean SED (i.e., the wavelength at which 
the EIR emission peaks shifts to shorter wavelengths); if 
it is added with a negative coefficient, the resulting SED 
is cooler than the mean. Thus, the coefficient of PCI can 
be considered a proxy for the effective dust temperature 
of the SED. In contrast, the second PC (PC2) affects the 
width of the SED: if it is added with a positive coefficient 
to the mean EIR SED, it tends to broaden the SED by 
increasing the power in the wings and removing power 
from the center. Conversely, it makes the EIR SED peak 
narrower if added with a negative coefficient. 

The coefficients of the PCs for a galaxy determine how 
its EIR SED differs from the mean SED. Eigure 2 shows 
how the coefficients of PCI and PC2 (we refer to them 
as Cl and C2, respectively) for each galaxy in the en¬ 
tire sample depend on four different global physical pa¬ 
rameters (SFR, AGN luminosity, IR luminosity and dust 
mass) that we expect to affect the FIR SEDs. We fo¬ 
cus on these specific parameters because they are simple 
global parameters that are clearly important for the ra¬ 
diative transfer, and they are typically available in SAMs. 
The first two characterize the radiation that can poten¬ 
tially heat the dust. The IR luminosity tells us how much 
radiation is absorbed by dust; of course, this quantity 
depends on both the bolometric luminosity of the stars 
and AGN and what fraction of the intrinsic luminosity 
is absorbed. The dust mass characterizes the radiation 
sinks. As detailed below, we have also explored using 
other individual or pairs of parameters to predict the co¬ 
efficient values, but we do not show them in this plot 
because we do not focus on them in the bulk of the anal¬ 
ysis below. 

Eigure 2 indicates that Cl correlates with IR luminos¬ 
ity (with a Pearson correlation coefficient of r = 0.71), 

Throughout this work, we use Ljr calculated by integrat¬ 
ing the SEDs over the wavelength range of 8 — 1000 /rm. When 
dust self-absorption is negligible, this quantity is independent of 
viewing angle and identical to the luminosity absorbed by dust. 
However, when dust self-absorption is non-negligible, Li^ can de¬ 
pend on viewing angle, whereas the luminosity absorbed by dust 
is an intrinsic property of the simulated galaxy that does not de¬ 
pend on viewing angle (see Hayward et al. 2011 for more details). 
We opt to use rather than the absorbed luminosity because 
(1) the former can be inferred from observations without recourse 
to radiative transfer modeling, (2) the average of Lir taken over 
a sufficient number of viewing angles must equal the absorbed lu¬ 
minosity, and (3) dust self-absorption does not lead to significant 
variation in Lir with viewing angle for the bulk of the simulated 
galaxies. 


































6 


M. Safarzadeh et al. 



Figure 3. Density plot of Cl versus C2. The color bar indicates 
the logarithm of the number of SEDs in the bin. Recall that higher 
values of Cl correspond to hotter SEDs, and higher values of C2 
correspond to broader SEDs. Cl and C2 are uncorrelated. 

SFR (r = 0.60) and AGN luminosity (r = 0.73), al¬ 
though there is a large scatter. This result indicates 
that as the IR luminosity, SFR or AGN luminosity is in¬ 
creased, the SED peak shifts toward shorter wavelengths, 
which is to be expected because of the known correla¬ 
tion between effective dust temperature and IR luminos¬ 
ity. The correlations between C2 and the four properties 
shown in Figure 2 are all weak negative correlations (all 
have —0.5 ^ < 0)- This result indicates that there is a 

weak tendency for higher-luminosity or/and higher-dust 
mass galaxies to have narrower FIR SEDs. 

Other than correlations between the PCs and the 
global parameters, it is also interesting to consider 
whether there is a correlation between the PC coeffi¬ 
cients. Although PCI and PC2 are orthogonal by con¬ 
struction, their coefficients may still be correlated. In 
our case, a correlation between Cl and C2 would indi¬ 
cate a relation between the effective temperature of the 
SED and its broadness. To see if there is such correla¬ 
tion, we plot Cl versus C2 for all the SEDs in our sample 
in Eigure 3. We find that Cl and C2 are uncorrelated. 

Because the scatter in the correlations between the co¬ 
efficients and global galaxy parameters shown in Figure 2 
is large, it is worth considering whether we can better 
predict the coefficients using two parameters simulta¬ 
neously. In particular, thermal equilibrium considera¬ 
tions suggest that combining dust mass and a luminosity- 
related parameter (e.g. IR luminosity or SFR) could be 
promising.Thus, in Figure 4, we show the results of 
predicting Cl (top) and C2 (bottom) in four different 
manners: 1) using logLiR (e.g.. Cl = AlogLi^ + B), 
2) using log SFR (e.g.. Cl = A log SFR +B), 3) using 
Lir and Mdust (e.g., Cl = AlogLiR -f RlogMdust + C) 
and 4) using log SFR together with logMdust (e.g., Cl 

We have not used the AGN luminosity as a predictor despite 
the correlation evident in Figure 2 for the following reason: for 
most of the simulated galaxies, the AGN contributes less than 10 
percent of the bolometric luminosity, and the FIR luminosity is not 
AGN-dominated (a detailed analysis of the contribution of AGN 
to FIR emission will be presented in Hayward et ah, in prep.) 
Gonsequently, the correlation between the AGN luminosity and 
G1 does not indicate causation; rather, it arises because in the 
simulations that we analyze, the black hole accretion rate and SFR 
are correlated. 



Predicted C2 - TrueC2 


Figure 4. The top panel shows the difference between the pre¬ 
dicted and true values of G1 (the coefficient of the first PC). Cl 
is estimated using five different estimators that depend on either 
one or two physical parameters; see the figure legend for the pa¬ 
rameters used. In all cases, the logarithm of the parameter is used. 
Using both IR luminosity and dust mass helps reduce the bias in 
the predicted Cl values compared with the other estimators. In the 
bottom panel, the same is shown for C2. Again, the combination 
of IR luminosity and dust mass is superior to the others, although 
the difference in predictive power is less than for Cl. 

= A log SFR -l-R logMdust + C)}^ It is clear that the 
combination of log Lir and log Mdust results is best able 
to predict Cl, and this combination is superior to log Lir 
alone (i.e., incorporating the dust mass reduces the error 
in the prediction). For C2, the combination of log Lir 
and logMdust is again superior to the others, but the 
differences in predictive power are less significant than 
for Cl. Our best-fitting relations for Cl and C2 are the 
following: 

Cl = 0,521„g - 0.471og (1^) - 1.88 (5) 

C2 = -0,0141og(ig)-0,101„g(^) (6) 

-fl.05. 

We note that the coefficients for log Lir and log Mdust in 
Equation (5) have similar magnitudes but opposite signs. 
This suggests that the value of Cl depends on the ratio 
Lir/M dust- We shall discuss this in detail in Section 8.1. 

We have explored using various other combinations of param- 
eters and found that these combinations had the best predictive 
power. 
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C2 depends only weakly on dust mass and is effectively 
independent of Fir. Thus, there is a mild tendency for 
sources with higher dust masses to have narrower SED 
peaks. Because of the weakness of the dependence, we 
will not interpret it further. 

We now investigate how well we can predict the SEDs 
by predicting the PC coefficients using Equations (5) and 
(fj) and using Equation (3) to predict the SED. Figure 
5 compares the predicted and true SEDs for 9 randomly 
chosen SEDs from our entire sample. The true SED is 
shown in blue. The red (green) line is the predicted SED 
obtained by adding the first PC (first and second PC) 
to the mean SED with the coefficient value(s) predicted 
based on the galaxy’s Lir and Mdust values (using Equa¬ 
tions 5 and 6) . The red (green) number in each panel in¬ 
dicates the reduced value (Equation 4) obtained when 
only PCI (both PCI and PC2) is (are) used to predict 
the SED. From these examples, we see that the FIR SEDs 
are generally predicted very well near their peaks. Gen¬ 
erally, when > 1, the reason is that the predicted and 
actual FIR SEDs differ at long wavelengths. Finally, us¬ 
ing PC2 improves the prediction (in particular, at long 
wavelengths) in only some cases, and it can actually make 
the prediction worse; this is true even if we use the true 
values of Cl and C2 rather than those predicted using 
Equations (5) and (6). We shall discuss these points in 
more detail below. 

Eigure 6 shows the distribution of Xr (Equation 4) for 
all SEDs in the full sample when the SEDs are predicted 
in this manner (the green line). We also show the Xr 
distributions obtained when using only PCI and predict¬ 
ing its coefficient. Cl, using either Lir alone (the blue 
dashed line) or Lir and M^ust (the red solid line). A com¬ 
parison of the solid red and blue dashed lines indicates 
that incorporating the dust mass results in significantly 
more accurate SED predictions compared with using Lir 
alone. The median Xr value is 0.56 (0.82) when both Lir 
and Mdust are (only Lir is) used to predict Cl. This re¬ 
sult indicates that the SEDs should be parameterized in 
terms of Lir and Mdust, not just Lir; we shall discuss 
this in detail below. The similarity of the green dashed 
and solid red lines indicates that using PC2 in addition to 
PCI does not yield significantly better SED predictions 
in terms of Xrl in this case, the median Xr is 0.54, only 
0.02 less than when PCI alone is used (and both Lir and 
Mdust are used to predict Cl). This is consistent with 
our observation from Eigure 5 that adding PC2 some¬ 
times leads to a better prediction at long wavelengths 
but sometimes causes the prediction to be worse. 

Eigure 6 indicates that for a subset of galaxies, the SED 
prediction fails catastrophically {Xr » 10)- explor¬ 
ing the locations of these galaxies in different parame¬ 
ter spaces, we determined that the catastrophic failures 
tend to have a high AGN contribution to the bolomet- 
ric luminosity or/and Lir > 10^^'® Lq, as indicated by 
Figure 7. This figure shows each galaxy in the plane of 
Lir and Lagn- The colors of the points indicate the 
X^ value of the predicted SED. The top panel shows the 
results for the prediction based on only PCI, whereas 
the bottom panel shows the results for the prediction 
based on using both PCI and PC2. Thus, the red cir¬ 
cles represent galaxies for which the SED prediction fails 
catastrophically. It is clear that in both cases, the galax¬ 


ies for which the SEDs are predicted least well tend to 
be galaxies that have high AGN luminosities given their 
Lir values (i.e., high AGN luminosity fractions) or/and 
have Lir > Lq. We analyzed the galaxies with 

high AGN fractions separately, but we were unable to 
determine global parameters that we could use to predict 
their SEDs well. We speculate regarding the reasons for 
our inability to predict the FIR SEDs of such galaxies 
in Section 8.4. We retain the high-AGN galaxies in our 
subsequent analysis, but the results are not significantly 
changed if we exclude them because only a small fraction 
of the simulated galaxies have Lagn >0.1. 

5. IMPACT OF DUST MASS ON THE SEDS 

EIR SEDs are often parameterized using templates 
that depend on the IR luminosity alone (e.g.. Chary & 
Elbaz 2001; Rieke et al. 2009; Lee et al. 2013; Symeonidis 
et al. 2013). At fixed redshift, the effective dust temper¬ 
ature is observed to increase (i.e., the peak of the EIR 
SED shifts to shorter wavelengths) with increasing Lir 
(C asey et al. 2014, and references therein). As discussed 
above. Figure 4 indicates that using both the IR lumi¬ 
nosity and dust mass increases our ability to predict the 
Cl coefficients for galaxies compared with using the IR 
luminosity alone as a predictor. Because we know that 
Cl basically makes an SED cooler or warmer with re¬ 
spect to the mean SED (higher Cl values tend to make 
the SEDs hotter), it is instructive to examine how Cl 
behaves on the Lir and Mdust plane. This is shown in 
Figure 8. One can see that at fixed Lir, higher values of 
Mdust correspond to lower values of Cl, i.e., cooler SEDs. 
At fixed Mdust, higher values of Lir lead to higher values 
of Cl, i.e., hotter SEDs. 

To more explicitly demonstrate the effect of the dust 
mass on the SEDs, we have re-run a subset of the dust 
radiative transfer calculations with artihcially lower and 
higher dust masses by changing the default dust-to-metal 
density ratio of 0.4 to 0.2 and 0.8, respectively. Con¬ 
sequently, the relative distribution of the dust remains 
the same, but the overall normalization and thus total 
dust masses are half or twice those of the standard run. 
Representative results for two snapshots are presented 
in Eigure 9; each panel shows how the FIR SED for a 
single time and viewing angle varies as the dust-to-metal 
density ratio is varied. The legends specify the assumed 
dust-to-metal ratio and values of Lir and Mdust for each 
of the SEDs. 

The results qualitatively agree with our expectations 
based on the above analysis. The top panel corresponds 
to the pre-coalescence phase of the most-massive major 
merger simulation from the z = 0 dataset. As the dust- 
to-metal ratio is increased, Lir increases, and the peak 
of the SED shifts slightly to longer wavelengths (i.e., the 
SED becomes colder). In this case, a substantial fraction 
of the luminosity is not absorbed when the dust-to-metal 
ratio is 0.2. Consequently, increasing the dust-to-metal 
ratio leads to higher optical depths and thus a larger 
fraction of light being absorbed and re-emitted in the 
IR. 

The SEDs shown in the bottom panel correspond to the 
M3M3e merger simulation near the peak of the starburst 
induced at final coalescence. Because the luminosity is 
dominated by a central dust-enshrouded starburst, the 
simulated galaxy is already opaque to effectively all of 
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Figure 5. Comparisons of the predicted and true SEDs for 9 SEDs randomly chosen from the entire sample. The blue line is the actual 
SED. The red (green) line is the predicted SED obtained by adding the first PC (first and second PC) to the mean SED with the coefficient 
value(s) predicted using Ljpr and Mjust. The numbers indicate the reduced value (Equation 4); a signal-to-noise ratio of 5 in each band 
was arbitrarily assumed. Note that using the second PC leads to a better prediction in only some of the cases. 


the radiation when the dust-to-metal ratio is 0.2. Conse¬ 
quently, Lir does not increase as the dust-to-metal ratio 
(and thus dust mass) is increased. Instead, it actually 
decreases by 0.1 dex; this occurs because we are consider¬ 
ing the Lir associated with a single viewing angle. When 
there is a non-negligible optical depth in the FIR, as can 
be the case for (U)LIRGs (see the discussion in Hayward 
et al. 2012), Lir can depend on the viewing angle. The 
increased viewing-angle dependence as the dust-to-metal 
ratio is increased explains the aforementioned decrease 
because as the dust-to-metal ratio is increased, more of 
the short-wavelength IR emission is removed from this 
particular line of sight. The absorbed luminosity, which 
is independent of viewing angle, is effectively identical in 
this case for all dust-to-metal ratios. 

Because the Lir values are almost identical for the 
three SEDs shown in the bottom panel, they provide 
a clean test of how the SED varies with dust mass for 
fixed IR luminosity. We see that as expected based on 
the trend shown in Eigure 8, the FIR SED systemati¬ 
cally shifts to longer wavelengths as the dust-to-metal 
ratio, and thus total dust mass, is increased; all other 
properties of the galaxy are kept fixed. 


The reason why at fixed dust mass, increasing the 
IR luminosity increases the temperature of the SED is 
that there are more photons available to heat the same 
amount of dust; consequently, thermal equilibrium dic¬ 
tates that the dust temperature must increase. The 
reason why at fixed IR luminosity, increasing the dust 
mass tends to make the SED cooler is that the luminos¬ 
ity is distributed over a greater mass of dust (because 
of dust self-absorption). Thus, thermal equilibrium im¬ 
plies that the dust temperature will decrease. Eor an 
isothermal modified blackbody, the temperature scales 
as T oc (TiR/Mdust)^^^^^^\ where /3 is the power-law in¬ 
dex of the dust opacity curve in the FIR (e.g. Hayward 
et al. 2011; Lanz et al. 2014). Although the simulated 
SEDs are not quantitatively well-described by this sim¬ 
ple model (Hayward et al. 2011, 2012; Lanz et al. 2014), 
it provides physical motivation for the claim that the ef¬ 
fective dust temperature increases (decreases) when Ljr 
(M dust) is increased and Mdust (Air) is kept fixed.Put 
another way, the mean intensity ‘seen’ by a dust grain 

For observational values of Lir and Mdust inferred from fitting 
galaxy SEDs with an isothermal modified blackbody, this relation 
is obeyed by construction. However, this is not the case for the 
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Figure 6. This figure shows how well the SEDs of the simulated 
galaxies can be predicted, as characterized by the Xr values, using 
the PCs. The red line corresponds to when only PCI is used; its 
coefficient is predicted using and M^ust- The dashed blue line 
shows the result when only PCI is used and its coefficient is pre¬ 
dicted using only The dot-dashed green line shows the result 

when both PCI and PC2 are used; their coefficients are predicted 
using both Lir and M^ust- When Cl is predicted using both Lir 
and Mdust? the Xr distribution shifts significantly to the left com¬ 
pared with when only Ljr is used; the median Xr value is 0.56 
(0.82) when both Lir and M^ust (only Tir is) used to predict 
Cl. Thus, incorporating the dust mass yields considerably better 
SED predictions. The similarity of the red solid and green dot- 
dashed lines indicates that using PC2 does not lead to significantly 
better Xr values, and the median value (0.54) is only slightly less 
than when only PCI is used (assuming that Cl is predicted using 
both Lir and Mdust)- 


is proportional to Ein/Mdust (see Draine & Li 2007 and 
Draine et al. 2007 for detailed discussions). Thus, as Lir 
is increased or Mdust is decreased, the mean intensity of 
light absorbed by the dust and thus typical grain tem¬ 
perature increases, and vice versa. 

6. IMPACT OF GALAXY SIZES ON THE SED SHAPE 

In addition to the total absorbed luminosity and dust 
mass, both of which must affect the SED shape because 
of thermal equilibrium, the geometry of radiation sources 
and dust can influence the SED. The surface densities of 
various components (e.g., all stars, young stars, gas, and 
dust) are global (observable) parameters that crudely 
characterize the global geometry of a galaxy and are often 
used in observational studies to interpret the evolution in 
the EIR SED shapes of galaxies (e.g., Elbaz et al. 2011; 
Rujopakarn et al. 2013). Thus, it is worthwhile to in¬ 
vestigate whether the PCA coefficients can be predicted 
better by incorporating information regarding the sizes. 

We calculated the 3D baryonic half-mass radii (re) and 
used these to approximate various average volume densi¬ 
ties by dividing integrated quantities - such as the SFR, 
Lir, and stellar, gas, and dust masses - by the half-mass 
radii cubed.We then investigated whether various 
combinations of volume densities could be used to predict 
the values of Cl and C2 better than using our standard 


simulations, in which the SED shape, and thus effective dust tem¬ 
perature, can in principle depend not only on the total luminosity 
heating the dust and the dust mass but also other factors, includ¬ 
ing the spatial distribution of the dust and sources and the dust 
composition. 

We opted to use volume densities rather than surface densities 
because the former are more relevant for the radiation transfer 
(albeit more difficult to infer from observations). 
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Figure 7. Shows the reduced x^ values for the predicted SEDs. 
In the top panel, the SEDs are predicted based on using only the 
first PC component, with its coefficient predicted using both IR 
luminosity and dust mass. In the bottom panel, the SEDs are 
predicted based on using both PCI and PC2; for both PCs, the 
coefficients are predicted using Ljr and dust mass, but the rela¬ 
tions used differ (see Equations 5 and 6). The colors of the circles 
indicate the value of the logarithm of the reduced x^ • The galaxies 
for which the SEDs are predicted least well tend to be galaxies 
that have high AGN luminosities given their Lir values or/and 
have Lir > L©. 


parameterizations, CI(LiR,Mdust) and C2 (Lir,M dust)- 

Figure 10 shows the difference between the true values 
of the PCA coefficients and those predicted using various 
estimators, as indicated in the legends (for clarity, we 
show only a subset of the combinations that we tried; 
the others fared comparably or worse). The top panel 
shows the results for Cl, and the bottom shows those for 
C2. To relate the coefficients’ values to a single galaxy 
parameter, we used models of the form Cl = A logX+B, 
where X is the galaxy parameter indicated in the figure 
legend and A and B are fitting coefficients. Similarly, 
when two parameters were used, we employed models of 
the form Cl = A log X + B log Y + C. 

Examination of Eigure 10 reveals that incorporating 
the size of the system when predicting the coefficients 
yields yields at best a very marginal improvement in the 
predictive power compared with using the total IR lumi¬ 
nosity and dust mass values alone. Thus, we have not 
used the sizes to predict the SEDs. This lack of improve¬ 
ment from incorporating the galaxy sizes suggests that, 
at least for the simulated galaxies, the overall spatial ex¬ 
tent does not play a significant role in determining the 
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Figure 8. Shows the value of the Cl coefficient (see the colorbar) 
as a function of Liji(on the y-axis) and M^ust (ori the x-axis). The 
interpolation between data points is performed using the nearest- 
neighbor method, which can result in some artifacts, but the trend 
that we see is robust. At fixed Lir, higher Mdust results in lower 
Cl values (i.e., lower effective dust temperature). At fixed Mdust? 
higher Ljr results in higher Cl values (i.e., higher effective dust 
temperature). 

shape of the FIR SED. We discuss this perhaps surprising 
result in detail in Section 8.3. 

7. TWO-PARAMETER FIR SED TEMPLATES 

Motivated by the results we have presented in the 
previous section, we now introduce a set of templates 
that is a function of both Lir and MdustMotivated 
by the PCA results, we separated the simulations’ FIR 
SEDs into (Fir, Mdust) bins and calculated the median 
SED in each bin. These are shown in Figure 11. Each 
panel shows the median FIR SED for each bin as a solid 
line, and the 16-84th percentile range is indicated by 
the shaded area. The dashed lines represent the SEDs 
predicted by adding PCI to the mean SED with the coef¬ 
ficient value predicted by inputting the median Ljr and 
Mdust values for each bin into Equation 5. The dot- 
dashed lines indicate the SEDs predicted when both PCI 
and PC2 are used. The SEDs have been normalized by 
dividing by Fir and are colored according to the wave¬ 
length at which the FIR peak is located, with redder col¬ 
ors indicating longer wavelengths. The trends revealed 
in the above analysis are apparent from the templates: 
(1) at fixed M^, the SEDs become warmer as Fir is in¬ 
creased (i.e., from left to right). (2) At fixed Fir, the 
SEDs become cooler when the Mdust is increased (i.e., 
from bottom to top). 

Figure 12 shows the ratios of the predicted SEDs to the 
true SEDs, F^/Fa, both when only PCI is used (solid 
lines) and when both PCI and PC2 are used (dashed 
lines). These figures demonstrate that the median SEDs 
in each bin are best predicted near their peaks. Short- 
ward of the SED peaks, the median SEDs of most bins 
can be predicted to within a few tens of percent. How¬ 
ever, for some bins, the prediction is worse; in one 
(the bin with median values Mdust = 10®'^ Mq and 
Fir = 10®'^ Lq), the SED is underpredicted by more 
than an order of magnitude at A ~ 35 /rm. At wave¬ 
lengths > 150 /rm, the SEDs are sometimes predicted 

The templates are available at the following URL: https: 
//www.cfa.harvard.edu/~chayward/fir_sed_templates.html. 




Figure 9. Shows the effect of changing the dust-to-metal ratio, 
and thus dust mass, on the SED. Each panel shows the SEDs of a 
simulated galaxy at a single time and viewed from a fixed viewing 
angle for three different dust-to-metal ratios, as specified in the 
legend; all other properties of the galaxies are unchanged. The 
Lir and Afdust values for each SED are also shown in the legend. 
In the top panel, the galaxy is not fully opaque when the dust- 
to-metal ratio is 0.2. Consequently, as the dust-to-metal ratio is 
increased, Lir increases by a factor of ~ 2, and the SED shifts 
only slightly to the right. In the bottom panel, Ljr decreases by 
0.1 dex as the dust-to-metal ratio is increased; this is a consequence 
of the viewing-angle dependence of Ljr (see the text for details). 
The SED systematically shifts to longer wavelengths as the dust- 
to-metal ratio is increased, which is consistent with our expectation 
based on the relationship shown in Figure 8 and considerations of 
thermal equilibrium. 

to within a few tens of percent. However, in many bins, 
the SEDs are underpredicted considerably. Using PC2 in 
addition to PCI sometimes makes the prediction better 
(i.e., the dashed line is closer to 1 than is the solid line). 

The goal of PCA is to explain the variance in a dataset. 
Thus, it is unsurprising that the peaks of the SEDs are 
predicted best, because this is the region that dominates 
the variance in the dataset. Because the luminosity den¬ 
sity at long wavelengths is ~ 1 — 2 orders of magnitude 
less than that at the SED peak, the long-wavelength re¬ 
gions of the SEDs contribute relatively little to the vari¬ 
ance. Consequently, it would be necessary to use higher- 
order PCs to predict the long wavelength emission well. 

8. DISCUSSION 

8.1. The simplicity of the FIR SEDs of galaxies 

We have found that the FIR SEDs of our simulated 
galaxies can be well predicted based on the galaxies’ 
IR luminosities and dust masses alone. Our results ex- 
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Figure 10. The difference between the predicted and true values 
of the coefficients for various estimators (similar to Figure 6). In 
this case, some of the estimators used include the galaxy sizes. The 
estimator used in each case is labeled in the legend, and in all cases, 
we use the logarithm of the parameter. For both Cl {top) and C2 
{bottom)^ incorporating the size does not significantly increase the 
predictive power. 

tend those of Hayward et al. (2011), who demonstrated 
that the observed-frame submm flux densities of simu¬ 
lated z ~ 2 SMGs could be well predicted using these 
two parameters. Moreover, as noted in Section 4, Cl 
effectively depends on logLm/Mdust because the best¬ 
fitting coefficients for (Lir and Mdust) have almost the 
same magnitude but opposite signs (0.52 and —0.47, re¬ 
spectively). Interestingly, in the simple isothermal mod¬ 
ified blackbody model, the dust temperature scales with 
Liu/Mdust-^^ In the more realistic case of a continuum 
of dust temperatures, the mean intensity of the radia¬ 
tion absorbed by the dust is proportional to FiR/Mdust 
(Draine & Li 2007). Thus, both simple models and our 
simulations suggest that the ratio Lm/Mdust is a key de¬ 
terminant of the FIR SED of a galaxy. 

These results indicate that the FIR SEDs of our sim¬ 
ulated galaxies are rather simple. The skeptical reader 

However, it is important to note that how the effective 
dust temperatures of the simulated galaxies’ SEDs depend on 
Tia/Mdust is not fully captured by this model (Lanz et al. 2014). 


may suggest that this simplicity is a consequence of the 
simplicity of our simulations. However, even in the sim¬ 
ulations, the FIR SEDs could in principle exhibit greater 
complexity because the simulated galaxies contain dust 
with a continuum of temperatures (which are set by the 
3D interstellar radiation field, 3D distribution of dust, 
and grain properties). Thus, the fact that Lm/Mdust en¬ 
codes so much of the variance in the simulated SEDs is 
somewhat surprising. 

Similar conclusions have been obtained based on ob¬ 
served FIR SEDs of galaxies. In particular, by htting 
observed galaxy SEDs using the model of Draine & Li 
(2007), Magdis et al. (2012) have argued that the redshift 
evolution of ‘main sequence’ galaxies’ SEDs is driven 
by redshift evolution in the Lm/Mdust ratio. For this 
reason, they suggest using SED templates for ‘main se¬ 
quence’ galaxies that depend on LiR/Mdust- 

Magnelli et al. (2014) studied how the effective temper¬ 
ature of galaxies depends on their position in the SFR- 
M* plane and redshift. They found that at all redshifts, 
the effective dust temperature smoothly increases with 
Lir, specific SFR, and the distance from the ‘main se¬ 
quence’ (i.e. excess SSFR relative to what one would 
expect for a ‘main sequence’ galaxy of the same mass), 
with the latter two correlations being more significant 
than the first. They interpret the dependence on the dis¬ 
tance from the ‘main sequence’ in terms of changes in the 
global star formation efficiency, SFR/Mgas. However, we 
note that for an approximately constant dust ratio, this 
quantity would serve as a proxy for LiR/Mdust- Thus, 
the results of Magnelli et al. (2014) are likely consistent 
with those of Magdis et al. (2012) and this work. 

8.2. Observational support for the importance of dust 

mass 

Observations indicate that at fixed Lir, z ~ 2—3 galax¬ 
ies have lower effective dust temperatures than do local 
galaxies (Casey et al. 2014 and references therein). Our 
results suggest that the observed trend could be a nat¬ 
ural consequence of high-redshift galaxies having more 
dust per unit IR luminosity. Thus, the evolution of the 
effective dust temperature-IR luminosity relation might 
be a consequence of evolution in the global properties of 
the ISM of galaxies rather than changes in the small-scale 
geometry of star-forming regions. 

The dust mass of a galaxy depends on the gas mass, 
gas-phase metallicity, and dust-to-metal ratio. As stars 
are formed, the ISM is enriched, which can increase the 
dust mass. However, star formation simultaneously re¬ 
duces the dust content of the ISM because the stars are 
formed from dust-enriched gas. Simple models that en¬ 
capsulate this competition between gas enrichment and 
consumption indicate that the maximum dust mass of a 
galaxy depends weakly on the gas fraction: it is max¬ 
imal when the gas fraction is 37 percent, but it varies 
by less than a factor of 3 for gas fractions in the range 
of 4-86 percent (Edmunds & Eales 1998). This result is 
for a closed box, but the limit also holds if outflows or 
unenriched inflows are allowed. Thus, although z ~ 2 —3 
ULIRCs are less metal-rich than local ULIRCs, they may 
still have higher dust masses. 

There is some observational evidence that supports our 
claim that the lower effective dust temperatures of z ~ 
2 — 3 ULIRCs are associated with higher dust masses 
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Figure 11. SEDs of our simulated galaxies binned according to IR luminosity and dust mass. In each panel, the median FIR SED in that 
bin is indicated by the solid line, and the shaded region represents the 16 - 84th percentile range of the FIR SEDs in that bin. The dashed 
(dot-dashed) line shows the SED predicted using PCI (PCI and PC2), with the coefficient(s) predicted using the median Liji and M^ust 
values for that bin. The SEDs have been normalized by dividing by Ljr. The color coding is based on where the peak of the FIR SED 
is located, with redder colors corresponding to longer peak wavelengths (i.e., colder effective dust temperatures). The logarithms of the 
median dust mass and IR luminosity (in solar units) are indicated in each panel, and the number of SEDs in each bin (N) is also shown. 
The dust mass increases in the upward direction, and Lir increases to the right. In most columns, the SEDs become redder from bottom 
to top; this visually illustrates the trend for increasing Af^ust to result in cooler SEDs when Lir is fixed. In a given row, the SEDs become 
bluer from left to right because at fixed M^usti the SEDs become hotter as Ljr is increased. The trend with Ljr is more apparent that 
the with Mdust because except for the bottom row, the Lir values in a given row span 2-3 orders of magnitude, whereas the M^ust values 
in a given column span only 1-2 orders of magnitude. Generally, the SEDs are predicted very well near their peaks. In some cases, the 
SEDs are underpredicted significantly at long wavelengths, although use of PC2 reduces the amount by which the SED is underpredicted 
compared with using only PCI. 
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Figure 12. The predicted FIR SED divided by the true FIR SED when PCI (PCI and PC2) is (are) used to predict the median SED 
in each bin is shown by the dashed (dot-dashed) lines. The solid lines correspond to a ratio of unity. The (TiR,Mdust) bins are the same 
as in Figure 11. This figure clearly shows that near their peaks, the median SEDs can be predicted well based on their Lir and M^ust 
values alone. At wavelengths significantly shorter or longer than the peak wavelength, the SEDs are predicted very well in some bins but 
are sometimes incorrect by greater than an order of magnitude. Incorporating PC2 makes the predictions more accurate in some cases, 
but it can also make them worse. 
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Figure 13. IR luminosity versus dust mass for two samples of IR- 
luminous galaxies, local LIRGs and ULIRGs from U et al. (2012) 
and 2 ~ 2—3 SMGs from da Cunha et al. (2015). The stars indicate 
the median dust mass for the subset of each sample within the 
shaded region, which is defined by 10^^'^ < Lir/ Lq < 1012-5. The 
median dust mass of the SMGs is an order of magnitude greater 
than that of the local (U)LIRGs. 

relative to local ULIRGs. Figure 13 shows Lir versus 
Mdust for a sample of local LIRGs and ULIRGs (U et al. 
2012) and a sample of z ~ 2 — 3 SMGs (da Gunha et al. 
2015). For Lipj ~ iq1i.5-12.5 local ULIRGs have 

a median dust mass of 10^-^^ Mq, whereas the z ~ 2 — 3 
SMGs have a median Mdust ~ 10®-®“^ M©. Thus, the 
SMGs have dust masses that are more than an order of 
magnitude greater than those of the z ~ 0 ULIRGs. We 
argue that the increased dust masses are the reason that 
the SMGs have cooler SEDs. This is consistent with the 
claim of Magdis et al. (2012), who argued that z ^ 0.5 — 2 
‘main sequence’ ULIRGs have cooler SEDs than local 
ULIRGs because the former have lower Lir/M dust ratios. 

We have used the Lir and Mdust data from da Gunha 
et al. (2015) because they were inferred from high- 
resolution ALMA data and Herschel data that were de- 
blended based on the ALMA data. Thus, their pho¬ 
tometry should be much less affected by blending than 
typical datasets from Herschel and other single-dish 
FIR/(sub)millimeter telescopes. This is highly desirable, 
because the fact that blending becomes more severe at 
longer wavelengths could cause the SEDs extracted from 
blended data to be colder than the true SEDs of indi¬ 
vidual sources. However, the significant caveat regard¬ 
ing our use of this dataset is that the SMG selection 
is biased toward colder effective dust temperatures, and 
thus the Liu/Mdust values of z ^ 2 SMGs likely do not 
represent the full range of Lm/Mdust values exhibited 
by z ~ 2 ULIRGs. To conclusively determine how the 
Lir/M dust ratios of z ~ 2 ULIRGs compare with those of 
local ULIRGs, high-resolution FIR and (sub)millimeter 
observations of a sufficiently large, unbiased (in terms of 
effective dust temperature) sample of z ~ 2 ULIRGs are 
required. 

8.3. The unimportance of galaxy sizes in determining 
the SED shape 

In Section 6 , we demonstrated that incorporating in¬ 
formation regarding the galaxy sizes did not significantly 


improve our ability to predict the FIR SEDs. This result 
may be surprising to some readers, given that it is often 
claimed that z ~ 2 galaxies have lower effective dust tem¬ 
peratures at fixed Lir because they are more extended 
(e.g., Elbaz et al. 2011; Rujopakarn et al. 2013). It is 
true that for a central source surrounded by dust (i.e., 
the ‘shell’ geometry of Misselt et al. 2001), increasing 
the spatial extent of the absorbing material will result 
in colder dust because the dust grains receive a more 
‘diluted’ radiation held. This geometry may be a rea¬ 
sonable approximation (especially if the shell is allowed 
to be clumpy) for individual HII regions and highly ob¬ 
scured AGN. Indeed, more compact HII regions exhibit 
hotter dust temperatures (e.g.. Groves et al. 2008). How¬ 
ever, the overall geometry of both real galaxies and our 
simulated galaxies is likely more similar to the ‘dusty’ 
geometry of Misselt et al. (2001), in which the stars and 
dust are mixed, because both the SFR density and dust 
density are correlated with the gas density (see also Jon- 
sson et al. 2006). In such a geometry, the temperature 
of the dust is insensitive to the size (assuming that the 
sizes of the stellar and dust distributions are scaled in 
the same manner; Misselt et al. 2001). 

8.4. The origin of catastrophic failures in the SED 
prediction 

Figure 7 indicates that most galaxies for which the 
SED prediction is a catastrophic failure (i.e., Xr » 10) 
have either high Lagn/Air or/and Lir > 10 ^^-® L 0 . 
There are a few potential reasons that the SEDs of such 
sources would prove to be especially difficult to predict. 
For the simulated galaxies in which the AGN contributes 
signihcantly to the bolometric luminosity, the AGN can 
heat host-galaxy dust and cause FIR emission (this will 
be discussed in detail in Hayward et al., in prep, and 
Roebuck et al., in prep.) Sources in which the AGN dom¬ 
inates the dust heating can differ from star formation- 
dominated sources in terms of the SED of the radiation 
absorbed by the dust. The geometry of such sources 
may also qualitatively differ: when the AGN dominates, 
the geometry is more similar to the ‘shell’ geometry of 
Misselt et al. (2001) than the mixed geometry, whereas 
the latter should better describe sources in which star 
formation dominates the dust heating, as argued in the 
previous subsection. Finally, in the AGN-dominated 
and most-IR luminous sources, dust self-absorption is 
likely more significant than in the other sources. Some 
or all of the above differences between star formation- 
powered and AGN-powered IR sources may explain why 
the SEDs of many of the sources with high Lagn/Lir 
or/and Lir > 10 ^^-® Lq cannot be predicted well. 

8.5. Implications for IR counts in hierarchical models 

Cosmological galaxy formation models have long strug¬ 
gled to correctly reproduce the observed IR and sub-mm 
counts without introducing fairly radical assumptions 
such as an extremely top-heavy stellar initial mass func¬ 
tion (Devriendt & Guiderdoni 2000; Baugh et al. 2005; 
Lacey et al. 2010; Dave et al. 2010; Somerville et al. 2012; 
Niemi et al. 2012, but cf. Hayward et al. 2013b; see also 
the discussion in Casey et al. 2014). Due to the infeasibil¬ 
ity of carrying out full 3D radiative transfer calculations 
on a cosmological hydrodynamic simulation (see 8 . 6 ), to 
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date most such calculations have relied on semi-analytic 
models combined with a simplified approach to comput¬ 
ing the FIR SEDs. For example, the models of Devriendt 
& Guiderdoni (2000) and Somerville et al. (2012) used 
empirical libraries of dust emission templates parame¬ 
terized only by Fir. These libraries clearly cannot cap¬ 
ture the observed redshift evolution of the relationship 
between effective dust temperature and IR luminosity. 
The SAMs presented by Somerville et al. (2012) and fur¬ 
ther investigated by Niemi et al. (2012) underpredicted 
IR number counts at wavelengths > 100 /xm, and the 
discrepancy became worse with increasing wavelength. 
Our work here suggests a straightforward and physically 
motivated way to improve the modeling of dust emis¬ 
sion in SAMs by using templates that depend on both 
Lir and Mdust- Our results are also encouraging for the 
use of SAMs to model dust emission for large samples 
of galaxies, as they suggest that the sizes and detailed 
geometries of galaxies (which are properties that SAMs 
cannot model accurately) are sub-dominant compared to 
the global parameters Ljr and Mdust- 

If indeed the dust is colder (at hxed Lir) in high- 
redshift galaxies, as indicated by observations, adopting 
two-parameter templates like the ones we have presented 
here will clearly work in the direction of alleviating the 
tension between the SAM predictions and observations. 
Accounting for the effects of blending will further reduce 
the discrepancy (Hayward et al. 2013a,b; Munoz Aran- 
cibia et al. 2015; Cowley et al. 2015). However, whether 
SAMs will predict a strong enough evolution of Mdust 
with Lir and redshift to reproduce the observed IR- 
submm counts, once observational effects such as blend¬ 
ing have been taken into account, remains to be seen. We 
plan to investigate this by incorporating two-parameter 
dust emission templates in the Somerville et al. (2008) 
SAM in a future work (Safarzadeh et al., in prep.). 

8.6. Limitations and future work 

The detailed simulation methodology that we have em¬ 
ployed in this work has the advantages of being well stud¬ 
ied, and in previous works, the specihc simulations used 
in this work and similar simulations have been demon¬ 
strated to reproduce the properties of a wide range of real 
galaxies, as discussed in Section 1. However, there are 
naturally some limitations. First, because of the manner 
in which our sample was constructed, the demograph¬ 
ics of the population are by no means representative 
of those of the real Universe. To achieve a cosmologi- 
cally representative population, it would be necessary to 
perform 3D radiative transfer on galaxies selected from 
a large-volume cosmological simulation. Unfortunately, 
such simulations typically have spatial resolution > 1 
kpc and thus do not resolve galaxy disk scaleheights, 
much less the internal structure of the ISM. Moreover, 
state-of-the-art cosmological simulations lack starbursts 
(i.e., there are significantly fewer outliers above the ‘star 
formation main sequence’ than observed; Sparre et al. 
2015), which are thought to power local ULIRGs and 
a non-negligible fraction of z ^ 2 — 3 ULIRGs (Hop¬ 
kins et al. 2010; Hayward et al. 2013a,b; Gowley et al. 
2015). Consequently, the utility of such simulations 
for investigating the FIR SEDs of galaxies remains lim¬ 
ited. Cosmological ‘zoom-in’ simulations can be used 
to achieve orders-of-magnitude better spatial resolution; 


thus, radiative transfer can be meaningfully applied to 
such simulations (e.g., Granato et al. 2015). However, 
the considerable computational expense of such simu¬ 
lations strongly constrains the subset of the parameter 
space that can be sampled. 

Even for idealized, comparatively simple simulations 
such as ours, the computational expense required to per¬ 
form the radiative transfer is significant. Consequently, 
the parameter space spanned by our simulation suite is 
not exhaustive, and the sampling of the parameter space 
is rather coarse. This limitation can be addressed in the 
future through the use of a (considerably) larger simu¬ 
lation suite or through performing radiative transfer on 
all resolved galaxies in a large-volume cosmological sim¬ 
ulation; however, we again stress that for the latter, the 
limited spatial resolution will continue to be a significant 
hurdle for the foreseeable future. We have no reason 
to expect that our qualitative conclusions would differ 
if we were to use a larger or/and cosmologically repre¬ 
sentative simulation suite, especially given the demon¬ 
strated agreement between our simulated galaxies’ and 
real galaxies’ SEDs. However, it is possible that the de¬ 
tails, such as the variation in the SED templates, are 
sensitive to the specific simulations used. 

Another significant limitation is that the hydrodynam- 
ical simulations do not resolve the detailed structure of 
the ISM, both because of the spatial resolution and the 
ISM model employed. As discussed in Section 2, we 
assume that the dust is uniformly distributed on sub¬ 
resolution scales. Using our current methods, it is pos¬ 
sibly to crudely characterize the uncertainty associated 
with the sub-resolution ISM structure by comparing two 
extremes, the default model and one in which the dust 
in the cold clouds implicit in the Springel & Hernquist 
(2003) ISM is completely ignored (i.e., the clumps have 
a volume filling factor of zero). We have performed such 
a comparison and found that the results do not quali¬ 
tatively differ. In fact, it is actually easier to predict 
the SEDs when the latter treatment is used. We spec¬ 
ulate that the reason for this result is that all optical 
depths are smaller (by construction); consequently, dust 
self-absorption is less significant. 

Moreover, we by no means claim that our treatments 
of stellar and AGN feedback are state-of-the-art, and 
various groups are now utilizing more sophisticated and 
likely more realistic feedback models (e.g., Agertz et al. 
2013; Agertz & Kravtsov 2015; Hopkins et al. 2014). 
However, to our knowledge, no UV-mm SEDs for such 
simulations have been presented in the literature; thus, 
the SEDs used in this work still represent the state-of- 
the-art for UV-mm SEDs computed by performing dust 
radiative transfer on hydrodynamical simulations. Com¬ 
putations of UV-mm SEDs of galaxies from the Feedback 
in Realistic Environments (FIRE) cosmological zoom-in 
simulations (Hopkins et al. 2014) using Sunrise are un¬ 
derway, but this is a significant undertaking in and of 
itself. Our method could be applied to these and other 
SED datasets in the future. 

Finally, we have demonstrated that PGA is a useful 
tool for identifying which parameters drive the variation 
in galaxy SEDs. However, the PGA results cannot be 
used to predict SEDs of galaxies outside of our parame¬ 
ter space because the mean SED and PCs depend on the 
dataset on which the PGA is performed. Moreover, be- 
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cause of the dependence on the mean SED and the fact 
that our simulation suite is not cosmologically represen¬ 
tative in terms of the galaxy demographics, the PCA 
results cannot even be used to predict the SEDs for sam¬ 
ples in which the parameter space is a subset of that 
spanned by our simulations but the distribution within 
the parameter space differs. Thus, the templates that 
we provide are a better tool for predicting FIR SEDs 
than are the PCA results. Moreover, other statistical 
methods, such as neural networks, may prove to be more 
useful than PCA for this purpose (e.g., Silva et al. 2012). 

9. CONCLUSIONS 

We performed PCA on a sample of FIR SEDs of simu¬ 
lated galaxies that we generated by performing dust ra¬ 
diative transfer on hydrodynamical simulations in post¬ 
processing. Our goal was to determine what drives the 
variation in galaxies’ FIR SEDs. Our main conclusions 
are the following: 

• The PCA indicated that only two PCs are suffi¬ 
cient to explain 97 percent of the variance in our 
SED sample. The first component characterizes the 
peak of the SED, whereas the second characterizes 
the breadth of the peak. 

• The coefficient of the first PC, Cl, is correlated 
with the IR luminosity, SFR, and AGN luminosity. 
This result indicates that the SEDs are hotter when 
the IR luminosity, SFR, or AGN luminosity are 
greater. 

• Incorporating dust mass increases our ability to 
predict the value of Cl and thus the FIR SEDs. 
At fixed IR luminosity, increased dust mass leads 
to lower Cl values and thus cooler SEDs. 

• The coefficient of the second PC, C2, is weakly anti¬ 
correlated with IR luminosity, SFR, AGN luminos¬ 
ity, and dust mass. It can also be predicted using 
Lir and Mdust, but the dependences on both quan¬ 
tities are weak. Using the second PC improves how 
well the SEDs can be predicted in some cases but 
makes the predictions worse in others. 

• Examination of the catastrophic failures to recon¬ 
struct SEDs revealed that the bulk of such SEDs 
correspond to simulated galaxies with high AGN 
fractions or/and Lir > 10^^ ® L©. For this sample, 
we were unable to predict the PC coefficients and 
thus SEDs well. 

• Incorporating galaxy sizes does not improve our 
ability to predict the SEDs. 

• The above conclusions suggest that the redshift 
evolution in effective dust temperature (i.e., at 
fixed Air, 0^2 — 3 galaxies exhibit lower effective 
dust temperatures compared with z ^ 0 galaxies) 
is not a consequence of higher-redshift ULIRGs be¬ 
ing more extended, as is often claimed. Instead, 
our work suggests that this difference is driven by 
z ~ 2 — 3 ULIRGs having higher dust masses at 
fixed Air (because of their higher gas fractions than 
local galaxies), as suggested by some observations. 


• Because of the importance of dust mass in deter¬ 
mining the FIR SED shape, a two parameter set 
of IR SED templates that depend on both Air 
and Mdust should be superior to those that depend 
on Air alone. We have generated such a set of 
templates based on our simulated SEDs and made 
them publicly available. They should be useful for 
fitting observed galaxy SEDs and predicting galaxy 
SEDs in unobserved wavelength regimes, and they 
can be used to predict IR SEDs of galaxies in cos¬ 
mological simulations and SAMs as long as the lu¬ 
minosity absorbed by dust and dust mass can be 
estimated. 
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